Quarkonium mass splittings in three-flavor lattice QCD 
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We report on calculations of the charmonium and bottomonium spectrum in lattice QCD. We 
use ensembles of gauge fields with three flavors of sea quarks, simulated with the asqtad improved 
action for staggered fermions. For the heavy quarks we employ the Fermilab interpretation of the 
clover action for Wilson fermions. These calculations provide a test of lattice QCD, including the 
theory of discretization errors for heavy quarks. We provide, therefore, a careful discussion of the 
results in light of the heavy-quark effective Lagrangian. By and large, we find that the computed 
results are in agreement with experiment, once parametric and discretization errors are taken into 
account. 
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I. INTRODUCTION 

Quarkonium plays an important role in the application 
of QCD to hadronic physics. Early calculations of char- 
monium based on potential models gave strong support 
to the interpretation of these states as bound states of 
a new heavy quark [l], Q • Although physically the char- 
monium state is analogous to positronium, its historical 
role as a model system of QCD proved to be analogous to 
that of the hydrogen atom in quantum mechanics. The 
charmonium spectrum provided a simple example of how 
QCD works, made even more compelling with the subse- 
quent observation of the bottomonium states. 

Although the analysis of the quarkonium spectrum 
based on potential models is a significant triumph of 
QCD, an ab initio calculation based on lattice QCD, 
an approach that can simultaneously deal with light 
quarks would be even more satisfying. However, in lattice 
QCD, the continuum limit requires that am approaches 
zero, where a is the lattice spacing and m is the mass 
of the state. Quarkonium states and their constituent 
quarks are so heavy, that it is often impractical to use 
so small a lattice spacing that am is small. Both lat- 
tice NRQCD [3,|| and the Fermilab action 5] have been 
developed to treat heavy quarks in lattice gauge theory. 
Thus, the successful calculation of the spectrum of char- 
monium and bottomonium becomes a significant test of 
these techniques. 

This paper describes the current state of the quarko- 
nium spectrum based on the Fermilab approach to heavy 
quarks, using gauge configurations provided by the MILC 
Collaboration [6( that incorporate the effects of three 
light quarks: up, down, and strange. We have been 
studying the quarkonium spectrum using this formal- 



ism for some time ML starting on ensembles with two 
flavors of sea quark [§[, and new ensembles of configu- 
rations have become available during the course of the 
project. In this paper, we report on results with four 
lattice spacings from a » 0.18 fm to « 0.09 fm. We do 
not yet consider this work a definitive calculation using 
our approach. Results from two finer lattice spacings 
should become available in the future. As we detail be- 
low, the tuning of the bare valence heavy-quark masses, 
via the heavy-light meson spectrum , is not yet precise 
enough to give satisfactory answers to all the questions 
that have arisen. In the future, we expect to have better 
control of the heavy-quark masses. Nevertheless, already 
we can successfully reproduce important features of the 
quarkonium spectrum, and we consider this another im- 
portant testbed in which to assess the errors that arise 
from our treatment of the heavy quarks. Knowledge of 
these errors is also important for calculations of proper- 
ties of heavy-light mesons, for example those pertaining 
to semileptonic decays [l(| , leptonic decays [ll[ , and B-B 
mixing [l|, as well as the heavy-light spectrum [{|. 

Prior lattice QCD work on quarkonium with differ- 
ent sea-quark content has been reviewed by Bali [l3l ]. 
More recently, Dudek and collaborators have used the 
quenched approximation to explore decays [l4j, radia- 
tive transitions [l5| . and the excited-state spectrum [r| 
of charmonium. Meinel has used lattice NRQCD to com- 
pute the bottomonium spectrum at one lattice spacing 
on four ensembles with 2+1 domain- wall sea quarks [l^ |. 
The HPQCD Collaboration has calculated the quarko- 
nium spectrum on many of the same MILC ensem- 
bles used in our work, using lattice NRQCD for bot- 
tomonium [TH, [Til , and using hi ghly improved staggered 
quarks (HISQ) for charmonium [20j. Lattice NRQCD is 
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not very accurate for charmonium, and HISQ requires 
very small lattice spacings for b quarks [2l|. An advan- 
tage of the Fermilab method is that it allows the same 
treatment for both the charmed and bottom quarks. Pre- 
dictions of the cb spectrum have also been made, for the 
pscudoscalar B c using NRQCD b quarks and the charmed 
quark propagators from this project [22j, and also for the 
vector B* using NRQCD b quarks and HISQ charmed 
quarks [23 |. 

The plan of this paper is as follows. In Sec. UH we de- 
scribe our methodology, explaining the actions used for 
gluons, light sea quarks, and the heavy valence quarks. 
We describe, in detail, how to use the heavy-quark effec- 
tive Lagrangian to understand heavy-quark discretiza- 
tion errors in quarkonium masses. We also discuss the 
construction of the hadronic correlators and how they 
are fit to determine meson masses. Several broad issues 
inform the uncertainties (statistical and systematic), and 
they are discussed in Sec. Mil In Sec. IIV1 we show our 
results for the splittings between various states or com- 
binations of states. Where possible, we have organized 
the presentation of the mass splittings around individual 
terms in the heavy-quark effective action, which clarifies 
the approach to the continuum limit at available lattice 
spacings. Section [V] contains our conclusions and sugges- 
tions on ways to improve on this calculation. 



II. METHODOLOGY 



line, such as IS 1 or 1 3 P. In particular, 

M(1S) = |(M„ +3M^), (2.1) 
M(I3P) = i (M Xc0 + 3M Xcl + 5M Xa2 ) , (2.2) 

using charmonium for illustration. For brevity we usually 
write IP for 1 3 P. These spin-averages are sensitive to 
the leading term in a nonrelativistic expansion. Note 
that the IP and l x Pi (also denoted h c and hb) levels are 
nearly the same in nature, which can be explained by the 
spin-spin interaction's short range — S(r) in the context 
of potential models [2|. 

Complementary to the spin-averaged masses are spin- 
splittings that hone in on spin-dependent corrections [25L 
|26[. Below, we examine the hyperfine splittings 

M{nS m ? s )=M JJ ^-M^, (2.3) 

using charmonium 15 notation on the right-hand side. 
For the P states two combinations are of interest 

M(nP spin _ orbit ) = i (5M Xc2 - 2M Xc0 ~ 3M Xcl ) , (2.4) 
M(nP tcnBOI ) = i (3M Xcl - M Xc2 - 2M Xc0 ) , (2.5) 

again using charmonium 15* notation on the right-hand 
side. M(nSnFs) and M(nP te nsor) are sensitive to spin- 
spin interactions, and M(nP sp i n _ or bit) to spin-orbit inter- 
actions. 



In this section, we collect several sets of information 
needed to understand the results that follow. Scction lll Al 
defines the notation for different quarkonium states and 
splittings. Then we provide details of the lattice gauge 
configurations that we have used in Sec. Ill Bl Next, 
in Sec. Ill C[ we review the Fermilab method, discussing 
in detail how it can be understood via an effective La- 
grangian. This discussion provides a link between the 
lattice fermion action and the computed mass splittings; 
an important theme in this paper is to scrutinize our 
numerical results according to these theoretical expecta- 
tions. Last, we explain how we form correlation func- 
tions in Sec. Ill Dl and how we fit them to obtain masses 
in Sec. UTEl 



A. Notation 

In this paper, we use two notations for hadrons and 
their masses, both the standard names from the Particle 
Data Group and the spectroscopic notation n 2S+1 Lj, 
where S, L, and J are the spin, orbital, and total angular 
momentum, respectively, of the nth radial excitation. As 
usual, L = 0, 1, 2, . . . are denoted S, P, D, .... 

It is often convenient to discuss spin-averaged masses 
(and mass splittings) , which we indicate with a horizontal 



B. Configuration details 

These calculations have been carried out on lattice 
gauge configurations provided by the MILC Collabora- 
tion [|| , listed in Table HI They were generated via the 
R algorithm [27j with the one-loop Symanzik-improved 
Liischer-Weisz gluon action [28| combined with 2 + 1 fla- 
vors of sea quarks simulated with the asqtad action 12911 . 
The n/-dependent part of the one- loop couplings [301 ] 
became available only after the ensembles were gener- 
ated. We have used ensembles at four lattice spacings: 
a « 0.18, 0.15, 0.12, and 0.09 fm (also called in the text 
"extra-coarse", "medium-coarse", "coarse" and "fine" 
ensembles, respectively). The first four columns of Ta- 
ble|I]list the parameters of these ensembles, including the 
masses of the sea quarks, denoting the pair as ami/am s . 
The lattice scale of each ensemble with different sea quark 
masses was kept approximately fixed using the length 
r\ (3lL l32l | from the static quark potential. The absolute 
scale from the T 2S'-1S' splitting was determined on most 
of our ensembles by the HPQCD Collaboration [TsL [T9I] . 
Details on the r\ determinations can be found in review 
of other work on the MILC ensembles [3311 . C ombining 
this determination with more recent work [34], [35j , leads 
us to take the range r\ = 0.318^q'qq7 fm in this paper. 
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TABLE I: Run parameters and configuration numbers for the ensembles used to study charmonium and bottomonium r\, J/ip, 
T, h, xo, and xi states with relativistic operators, and h and all xj states with nonrelativistic operators. The labels a and b 
are used to distinguish between the two runs with the same j3 = 6.76 but different ami/am 3 . 













relativistic 






nonrelativistic 




a (fm) 


/3 


ami /arris N% 


x N t 




N c f 

- 1 * cont 


Kb 


N b c 

A * cont 


K c 


- 1 * cont 


Kb 


N b c 

A * cont 


« 0.18 


6.503 


0.0492/0.082 16 a 


x 48 


0.120 


401 






0.120 


400 









6.485 


0.0328/0.082 


11 


11 


331 






11 


501 








6.467 


0.0164/0.082 


n 


11 


645 






11 


647 








6.458 


0.0082/0.082 


11 


11 


400 






11 


601 






w 0.15 


6.600 


0.0290/0.0484 16 a 


x 48 


— 








0.122 


580 


0.076 


595 




6.586 


0.0194/0.0484 




0.122 


631 


0.076 


631 


11 


580 




595 




D.OiZ 


n nnQ7 /n n/is/i 


11 


11 


631 


11 


(\ 3 1 


11 


629 


)) 


DOl 




6.566 


0.00484/0.0484 20 3 


x 48 










11 


601 


11 


600 


w 0.12 


6.81 


0.03/0.05 20^ 


x 64 


0.122 


549 


0.086 


549 












6.79 


0.02/0.05 


n 


11 


460 


11 


460 












6.76, a 


0.01/0.05 


11 


11 


593 


11 


539 












6.76,6 


0.007/0.05 


11 


?j 


403 














w 0.09 


7.11 


0.0124/0.031 28 a 


x 96 


0.127 


517 


0.0923 


517 


0.127 


518 


0.0923 


510 




7.09 


0.0062/0.031 


n 


11 


557 


11 


557 


11 


557 


11 


557 




7.08 


0.0031/0.031 40 3 


x 96 


n 


504 


11 


504 


11 


504 


)) 


504 



C. Heavy quark formulation 



In this work, the charmed and bottom quarks are sim- 
ulated with the Fermilab action Q 



n n 

- f - / n V^n j 



(2.6) 



where ?7 denotes the gluon field, and ip and ip denote 
the quark and antiquark fields. The clover definitions of 
the chromomagnetic and chromoelectric fields B and E 
are standard and given, for example, in Ref. p|. When 
C = r s = 1 and cb = ce = 0, S reduces to the Wilson 
action 36]; when £ = r s — 1 and cb = ce = csw> it 
reduces to the Sheikholeslami-Wohlert action [37|- The 
relation between the hopping parameter k and the bare 



mass is 



m a = --l 



3r,C 



(2.7) 



in four space-time dimensions. 

To motivate our choices of the input parameters k, £, 
r s , cb, and eg, let us review the nonrelativistic interpre- 
tation of Wilson fermions [5]. The pole energy of a single 



quark of spatial momentum p is 



E(p) = mi 



2m 2 



(2.8) 



where the quark rest mass m\ and the kinetic mass 771 2 
are defined at all orders of perturbation theory via the 
self energy [38| . At the tree level, 



mi = a 1 ln(l + toocz), 

J_ = 2C 2 | ar s C 

?ri2 mo(2 + moa) 1 + mga 



(2.9) 
(2.10) 



In general, mi 7^ r«2 unless mod <C 1; for charmed and 
bottom quarks on the ensembles listed in Table U one has 
mQ C a < 1, mofca > 1. One could tune £ so that m2 = mi, 
and we shall revisit that strategy below. 

Equation (12. 8p is the simplest example of a nonrel- 
ativistic interpretation of physical quantities computed 
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with the action in Eq. (|2.6p . This is justified, because in 
quarkonium the relative momentum of the heavy quarks 
is small compared with the heavy-quark mass. This is the 
basis of the phenomenological success of potential mod- 
els, which yield estimates of the relative velocity and, 
equivalently, internal momentum. For charmonium 

v w 0.55, p w 840 MeV, (2.11) 

and for bottomonium 

v ~ 0.31, p w 1475 MeV. (2.12) 

For both systems the typical kinetic energy is 450 MeV, 
as seen, for example, in the IP-IS* splitting. The ki- 
netic energies ^iri2V 2 are small on our lattices, and the 
momenta m-iv are marginally small (especially for bot- 
tomonium) . 

The nonrelativistic interpretation can be extended be- 
yond the tree level and to higher order in the nonrelativis- 



tic expansion using effective field theories. This has been 
pursued in detail emphasizing heavy-light hadrons (39l . 
S3, EH , and here we explain the ideas in the context of 
quarkonium. As in the Symanzik effective theory, one 
introduces a continuum effective Lagrangian, but here it 
is an effective Lagrangian valid for heavy quarks. With 
quarkonium, the appropriate power-counting rule for the 
effective Lagrangian is that of nonrelativistic QED [42j 
and nonrelativistic QCD (NRQCD) 0,i, S|. So one has 

S = -J2[ ^4q, (2-13) 

where s counts the powers of velocity. Here = means that 
the lattice gauge theory on the left-hand side, defined in 
our case by Eq. (|2.6p . is given an effective description 
by the right-hand side. The first several terms of the 
effective Lagrangian are 



(2) _ 
HQ — 



r (4) 
^HQ 



-ft (+) (D 4 + mi)/i (+) 



2m B 
ia ■ Bh(- 



2m 2 
U+\D 2 ) 2 h^ 
+ 8ml 

U-^ur ■ (D x E)h^ U~\D-E)h^ h^{D 2 ) 2 h^ 



2m 2 

Ha ■ (D x E)h^ 
8m| 



U+\D-E)hM 



8m|, 



+ ±a 3 W4 fr (+) A 4 ^ (+) 



2m b 



8m| 



8m|, 



8mf 



r 



(2.14) 



+ \cPw i h i - ) Dfh { -\ (2.15) 



where is a two-component field describing the quark, 
and /i^ - ) is a two-component field describing the anti- 
quark. The short-distance coefficients mi, , Trig , 
3 , and u>4 depend on the bare quark 



"E ' 



"E' ' 



masses, the bare gauge coupling, and all other couplings 
of the (improved) lattice action. The terms in scale 
with the heavy quark's velocity as v s , with the rules [4] 
D ~ ni2V, E ~ m^v 3 , and B ~ m\v . In particular, 
the nonrelativistic kinetic energy, D 2 /2m,2 



i ? 



IS 



an essential part of quarkonium dynamics, which is why 
mi appears with v in the power counting. The short- 
distance coefficients ttiZ , m^, 2 , etc., can be expanded 
in perturbation theory, with a s ~ v [4]. We have put 
the rest mass mift^'/i^ and temporal kinetic energy 
h^Dih^ into £hq, because by the equation of motion 
£>4 + mi ~ D 2 /2m2 ~ \m2V 2 . The next set of terms, 

£hq, are not written out, because they are numerous yet 
merely describe subleading contributions to the splittings 
examined below. 

One would like to adjust k, £, r s , cb, and ce so that 
the lattice gauge theory matches continuum QCD with 
controllable uncertainty. One would also like to reduce 
the number of input parameters as much as possible, to 
make the simulation easier to carry out. The coupling r s 



is redundant: any choice is allowed as long as the dou- 
bling problem is solved. We take 



= 1. 



(2.16) 



To derive tuning criteria for the others, one refers to 
the NRQCD description of continuum QCD, which takes 
the same form as Eqs. (|2.14[) and (|2.15|) , but with the 
following substitutions: 



mi h 


-> m, 


(2.17) 


m 2 h 


-> m, 


(2.18) 


1 


Z B 
-> , 


(2.19) 


m_e 


m 




1 

— *~ 


Ze 

~* ~~ 2' 


(2.20) 


m% 






1 


Ze> 


(2.21) 


2 h 


if' 


m E> 


m* 


1 

S H 


Zi 

_ 3' 


(2.22) 


m| 


m° 




W4 H 


+ o, 


(2.23) 



where the last is a consequence of Lorentz invariance, as 
is the exact equality of the rest and kinetic masses. The 



5 



matching factors Zi are unity at the tree level and have 
a perturbative expansion. To bring the lattice field the- 
ory in line with continuum QCD, one must then simply 
adjust the lattice couplings so that the lattice quantities 
on the left in (|2.17[) - (|2.23[) become, to some accuracy, 
the continuum quantities on the right. In principle, this 
matching could be carried out nonperturbatively [44] , al- 
though we do not pursue that strategy here. 

If one restricts one's attention to mass splittings and 
matrix elements, it is not necessary to adjust a coupling 
to tunc m\. The operators h^h^ are number oper- 
ators, commuting with everything else in the Hamilto- 
nian (39|. It is therefore acceptable to tolerate a large 
discretization error in the rest mass, and, consequently, 
one does not need to adjust £. We take 



C = i- 



(2.24) 



To obtain the correct dynamics, one must adjust k so 
that the rest of £^q is correctly tuned. In other words, 
one must identify the kinetic quark mass mi with the 
physical quark mass. 

The adjustment of cb stems from a concrete realization 
of (|2~T!)1) . At the tree level 



1 

m B 



2C 2 



ac B ( 

mo (2 + mo a) 1 + moa' 



(2.25) 



so to ensure tub — vni (as desired at the tree level where 
Zb = 1), one needs c B =r s . In practice, we take [recall- 
ing Eq. (H5) ] 



cb 



(2.26) 



to account for tadpole diagrams at higher orders in per- 
turbation theory [45j. On the coarse ensembles, we set 
Mo from the Landau link; on the other ensembles, we set 
it from the plaquette. 

In principle, the adjustment of ce should stem 
from (|2.20p . These simulations have been carried out, 
however, in concert with calculations of heavy-light 
masses [9fl. for which the adjustment of ce is a subleading 
effect [39|, |46[. Thus, we have taken 



Using formulae in Ref. [5l| . we can estimate the error 
stemming from \/m E , finding a tree-level mismatch of 



1 



1 



Ami 



(2 + moa)(l + mod) 4(1 + moa) 2 ' 

(2.28) 

where the right-hand side holds for £ = r s = Cb = Ce = 
1. At the tree level = m^, so the same error is made 
in the Darwin terms hS^'D ■ Eh^\ 



ce = cb- 



(2.27) 



An advantage of using Eqs. (UdD, ([2~Ti)l and (|2~TS)l to 

describe our lattice calculation is that it clarifies which 
parameters in S play a key role in various splittings de- 
fined in Sec. Ill Al The spin-averaged masses receive en- 
ergy (beyond 27711) from the balance between the kinetic 
energies h^D 2 h^ and the exchange of temporal glu- 
ons between h^A^h^ and h^A^h^'. As discussed 
above, they are sensitive to 7712, motivating the tuning of 
k (and the fixed choice for ^.) The hyperfine splittings 
M (nSnFs) arise from exchange of spatial gluons between 
h^icr-Bh^ and h^hcr-Bh^\ Hence they are propor- 
tional to l/m B and, drilling further back to S, sensitive 
to the coupling cb- The same line of dependency holds 
for the tensor splittings M(nP tcnsor ). Similarly, the spin- 
orbit part of the XcJ and Xbj levels arise from exchange 
of a temporal gluon between h^icr ■ (D x E)h^ and 
Aih^\ Hence they are proportional to l/m E and, 
referring back to S 1 , sensitive to eg. 

With the tree-level adjustment of cb, the hyperfine 
splittings should be expected to have errors of order 
a s mv 4 from radiative corrections to m^ 1 [relative er- 
ror: 0(a s ) ~ O(u)], and of order v e from the terms 

U^{D 2 ,i(T-B}h^ m.d§ Q [relative error: 0(v 2 )}. Sim- 
ilarly, with ce — cb, we expect leading errors of order 
a 2 m 3 v in the spin-orbit part of the x splittings [relative 
error: 0(m 2 a 2 )], as well as radiative corrections to m E 2 
[relative error: again 0(a s ) ~ 0(h)] . On the MILC en- 
sembles both relative errors are expected to be a few 
to several percent [5l| . and, perhaps counterintuitively, 
smaller for bottomonium than charmonium [5l| . 

The lattice action in Eq. (|2.6p does not contain pa- 
rameters to tune the two terms proportional to p 4 in 
Eq. (|2.15[) . The mismatches 



a 2 a 2 (l + 4m a) m a 4 



and 



8ml 8ml 2mo(2 + moa) 2 (l + moa) 4mo(2 + moa)(l + moa) 2 ' 8(l + moa) 

I 

fl 3 w _ 2a 2 a 3 ^ ative errors, 0(m 2 a 2 ), are again expected to be a few 

mo(2 + moa) 4(1 + moa) ~ to several percent, but in this case larger for bottomo- 

nium than for charmonium [51] . For plots of the a de- 
(given again for ( = r s = c B = c E = I) cause errors of pen dence of discretization effects caused by Eqs. (|2~2^1) . 
order a 2 m 3 u 4 in the spin-averaged splittings. The rel- 
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(l2~29l) . and (12301) . see Figs. 2 and 3 of Ref. [5l|. 

To tune k nonperturbatively, one adjusts it so that 
a hadron mass agrees with the experimentally measured 
value. Let us define Mi and M 2 for a hadron analogously 
to Eq. (I2.8P - 1 From the effective Lagrangian description 
for quarkonium, Eqs. (|2.13p - (|2.15[) . it follows that 



Mi = 2mi + Bi, 
M 2 = 2m 2 + B 2 , 



(2.31) 
(2.32) 



where the binding energy Bi is determined by terms of 
order v 2 and higher, but B 2 by terms of order v and 



higher 52]. In the splittings of rest masses, mi drops 
out, so we can obtain well-tuned results for Bi (and 
their differences) by adjusting k so that m% corresponds 
to a physical quark. That suggests tuning k so that, say, 
M 2 (IS) agrees with experiment. The spin average is use- 
ful, because it eliminates the leading effect of a mistuned 
chromomagnetic coupling cb- 

A better approach, still using a hadron's kinetic mass, 
is as follows. Reference [HJ analyzes the Breit equa- 
tion to show how higher-order potentials and the p 
terms generate B 2l tracing how the mismatches noted 
in Eqs. and ([23!!]) propagate to B 2 . This analy- 

sis reveals that the discretization error in B 2 is smaller 
for heavy-light hadrons than for quarkonium states. For 
heavy-light hadrons, the largest part of the kinetic bind- 
ing energy comes from the light quarks and gluons, and, 
since the light quark has mass ma <C 1, its contribution 
to the kinetic binding energy of the meson has only a 
small discretization error. To tune k for charmed and 

bottom quarks, it is therefore better to use heavy-light 

(*) (*) 

states, such as Ds and Bg , whose kinetic masses have 
the smallest statistical, discretization, and chiral extrap- 
olation errors. In fact, the leading discretization error, 
from the chromomagnetic energy, can again be removed 
by taking the spin-averaged mass of the pseudoscalar and 
vector mesons. 

It is sometimes thought that the tuning inaccuracy of 
the kinetic binding energy B 2 can be circumvented by 
adjusting £ so that (a hadron's) Mi = M 2 , and then 
fixing Mi to experiment. But any discretization error in 
B 2 is then propagated to £ and, hence, throughout the 
rest of the simulation. It is, therefore, just as clean to 
leave £ = 1 and tune M 2 to a target meson mass, as we 
have done here. 

At this stage, it may be helpful to compare and 
contrast the Fermilab approach 0, [5l| with lattice 
NRQCD [1, H. The construction of lattice NRQCD 
starts with the (dimensionally regulated and MS- 
renormalized) NRQCD effective Lagrangian for contin- 
uum QCD [i^, IHj], an d then discretizes it. This pro- 
cess can be repeated order-by-order in perturbation the- 



ory. In the Fermilab method, a version of the Wilson- 
Shcikholeslami-Wohlert lattice action is used, but the re- 
sults are interpreted with (dimensionally regulated and 
MS-renormalized) NRQCD with modified short-distance 
coefficients. This is possible because Wilson fermions 
possess heavy-quark symmetry, and the proposed im- 
provements preserve this feature. Then the parallel struc- 
ture of the NRQCD descriptions of QCD and lattice 
gauge theory are used to match the latter to the former. 
In both frameworks, the lattice action can be systemati- 
cally improved via the nonrelativistic expansion 0, l5lll. 

At a practical level, early spectrum calculations [47| 
use a lattice-NRQCD action [4{ that adjusts, at the tree 
level, the full v A Lagrangian and the spin-dependent v 6 
Lagrangian. 2 The p A terms are, thus, correctly normal- 
ized at the tree level, so the quarkonium and heavy-light 
kinetic mass tunings are comparably accurate. On the 
other hand, the Fermilab action has tree-level errors in 
the v 6 and even some of the v terms. The errors diminish 
monotonically as a is reduced, however. This is especially 
important for charmonium: here the nonrelativistic ex- 
pansion is not especially good, but it is needed only to 
organize the matching of the most important couplings 
in S, knowing that further errors, such as those described 
by jChq, are of the form (mv 2 a) 2 and smaller. 

In summary, the pattern of discretization effects leads 
us to tune k via kinetic masses corresponding to the IS* 
D s and B s mesons. The main spectroscopic results, pre- 
sented in Sec. IIV1 are for mass splittings, in which case 
the uncertainties are minimized by quoting differences of 
our computed rest masses. 



D. Correlator construction 

The meson correlator at a given spatial momentum p 
and time t is defined as 

C ah (p,t) =J2 e ~ iP ' X (0\O a (x,t)Ot(0,0)\0), (2.33) 

x 

where x is the spatial coordinate. The source and sink 
meson operators Ob and O a have the form 

O c (x, t)=J2 t)T(b c (x - y)V(y, *), (2.34) 



where T is a product of Dirac matrices appropriate for 
the meson spin structure, and <f) c (x — y) is a smearing 
function. Neglecting the disconnected piece, the meson 
correlator can be re-written with the quark propagators 



G(x, t; 0, 0) = / [dtp] [dip]ip(x, t)i/j(0, 0)e 



(2.35) 



In this paper, we use mi, rri2, ■ ■ . for quark masses, and Mi, 
M2, ■ ■ ■ for hadron masses. 



2 The HPQCD Collaboration's most recently published un- 
quenched calculations |l9( of the bottomonium spectrum with 
lattice NRQCD are obtained from an action without the spin- 
dependent v e corrections. 
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TABLE II: Nonrelativistic meson operators for the IP states. 
The smearing operator in spatial direction i is denoted by pi. 
The indices j and k are different from i and each other, and 
repeated indices on the last line are not summed over. 



Meson 




Irrep. 


Operator 




h 


'Pi 


Ti 


Pi, i = 


1,2,3 


Xo 


3 P 


Ai 


Ei=l a iP* 




Xi 


3 Pl 


Ti 


Oj x pk, i = 


1,2,3 


X2 


3 P2 


T-> 


",./<; + r>; /;.. % = 


1,2,3 


X2 


3 P2 


E% 


a jPi — VkPk, i — 


1,2 



TABLE III: Prior central values for the ground-state masses. 
The priors' widths are all fixed to 0.5. 



a (fm) 


k 


M qq a 


« 0.18 


0.120 


1.932386 


wO.15 


0.122 


1.841549 




0.076 


3.818718 


« 0.12 


0.122 


1.539279 




0.086 


3.187431 


« 0.09 


0.127 


1.184840 




0.0923 


2.818421 



with S from Eq. (|2.6|) . yielding 

Cd(p,t)=£e-*'x (2.36) 

a; 

Tr[G(O,0;a;,i)rG a6 (a;,i;O,0)rt] , 

where 

G a6 (x, i; 0, 0) = £ a (^ - y)G(y, i; *, 0)0* (*) (2.37) 

is the smeared quark propagator. 

For the P states, we use two types of quarkonium corre- 
lators, which we call "relativistic" and "nonrelativistic." 
In the relativistic case, all four spin components of the 
quark propagators were used to construct the two-point 
functions. We used point and smeared sources and sinks. 
The smearing functions 4>c(x) are IS and 2S wavefunc- 
tions of the QCD-motivated Richardson potential 
At the sink, spatial momentum p = 2ir(ni, n 2, n$)/L is 
given to the quarkonium state. We restrict the range of 
p such that J2 n l < 9- Using this approach, we com- 
puted correlation functions for the IS and 25 states for 
the pseudoscalar and the vector to study both the ki- 
netic and rest masses. For the IP states h, xo an d Xx we 
computed only the rest masses. 

In the nonrelativistic approach to constructing the 
two-point functions, the meson operators project onto 
two of the Dirac components of the quark fields. Ta- 
ble [IT] gives the explicit form of these operators. At the 
source and sink we smear the quark propagators with a 
P-type wavefunction <j) c (r) — faq(\r\)fj where </>is(|r|) 
is a Richardson IS wavefunction [531 ] and i = 1,2,3. At 
the origin we set C (O) = 0. The relativistic interpo- 
lating operators include extra lower Dirac components 
that increase the overlap with excited states. Therefore, 
one should expect that the overlap of the nonrelativistic 
meson operators with the IP ground states to be better 
than in the relativistic case. We used these nonrelativis- 
tic operators at p — for the h, XQi Xi> an d Xi states. 
In Sec. IIII Al we compare the results for the first three 
states with the corresponding results from relativistic op- 
erators. 

For both correlator constructions, we use several time- 
slice positions for the source vectors. In the case of the 
coarse (3 = 6.76, ami/am s — 0.005/0.05 ensemble and 



all medium-coarse ensembles, we use eight sources for 
the relativistic operators; in all other cases, we use four. 



E. Fitting methods 

To determine the mass spectrum, we fit our correlator 
data with a Bayesian procedure, taking priors guided by 
potential models [H, The priors, listed in Table [TTTT 
are the same for both relativistic and nonrelativistic cor- 
relators. To find the quarkonium masses from relativistic 
correlators, we use a delta function and a 15 smearing 
wavefunction as the source and sink. We fit simultane- 
ously two or three source-sink combinations for the zero- 
momentum states, including the ground state and up to 
two excited states. The minimum and maximum source- 
sink separation is varied, and the best fit is selected based 
on the confidence level and the size of the errors in the 
ground state and first excited state masses. After choos- 
ing the fit range, 250 bootstrap samples are generated to 
provide an error estimate. 

The fitting method in the case of nonrelativistic op- 
erators is similar except we use the same P-type wave- 
function, described above, for both source and sink. In 
this case, we use no more than a ground state plus one 
excited state in the fitting form. The quality of data in 
the nonrelativistic case is such that often a fit with just 
the ground state is enough, provided the fitting range is 
appropriately chosen. 



III. GENERAL RESULTS 

Before presenting results for mass splittings (in 
Sec. lIVp we discuss three general issues: a comparison 
of the statistical quality of relativistic and nonrelativis- 
tic operators (Sec. lIIl"A")) ; a numerical comparison of tun- 
ing k via M2 in heavy- light and quarkonium (Sec. IIII B[) : 
and a discussion of how uncertainties from tuning k are 
propagated to the mass splitting (Sec. IIII Cj) . 
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TABLE IV: Rest masses of the charmonium states r] c , J/ip, h c \c0, and xa calculated with relativistic operators. All masses 
in units of n = 0.318 fm. The star denotes masses that differ from their counterparts in Table IVl by more than 1.5a. 



a (fm) 


/3 






Vc^ 1 So) 




V(2"S'i) 




Xco(l J Po) 


Xd(l J A) 


» 0.18 


6.503 
6.485 
6.467 
6.458 


0.120 

5) 
5) 


3.2924(9) 
3.3071(14) 
3.3327(7) 
3.3481(13) 


4.24(6) 
4 42(4) 
4.39(27) 
4.47(6) 


3.4452(16) 
3.4581(18) 
3.4862(11) 
3.5004(16) 


4.35(7) 
4.48(3) 
4.45(11) 
4.43(10) 


4.185(17) 
4.214(26) 
4.213(29) 
* 4.217(18) 


★ 4.079(12) 
4.117(15) 
4.109(12) 
4.106(18) 


★ 4.052(89) 

★ 4.173(13) 
4.200(25) 

★ 4.181(19) 


» 0.15 


6.586 
6.572 


0.122 

If 


3.5688(8) 
3.5883(9) 


4.66(3) 
4.64(5) 


3.7317(13) 
3.7501(14) 


4.75(3) 
4.79(2) 


★ 4.476(8) 

★ 4.495(8) 


4.341(6) 
4.368(6) 


4.450(7) 
4.471(17) 


» 0.12 


6.81 
6.79 
6.76, a 
6.76,6 


0.122 

If 
J) 
53 


3.8721(11) 
3.8876(12) 
3.8824(9) 
3.9009(8) 


5.16(4) 
5.14(3) 
5.09(4) 
5.12(3) 


4.0594(18) 
4.0747(18) 
4.0677(15) 
4.0864(11) 


5.25(3) 
5.22(4) 
5.10(5) 
5.27(3) 


4.807(17) 
4.821(12) 
4.800(13) 
4.817(14) 


4.626(10) 
4.657(10) 
4.658(8) 
4.650(30) 


4.755(19) 
4.791(10) 
4.758(15) 
4.785(15) 


» 0.09 


7.11 
7.09 
7.08 


0.127 

5) 
5) 


4.2740(26) 
4.2885(15) 
4.2889(26) 


5.33(13) 

5.52(4) 

5.51(5) 


4.4460(22) 
4.4596(15) 
4.4613(33) 


5.55(6) 
5.66(4) 
5.65(7) 


5.159(29) 
★ 5.149(24) 
5.167(33) 


5.027(16) 

★ 4.986(15) 

★ 4.986(26) 


5.185(10) 
★ 5.123(19) 
5.133(48) 



TABLE V: Rest masses of the charmonium states h c , XcO, X^l, and x<=2 calculated with nonrelativistic operators. All masses 
in units of n = 0.318 fm. The star denotes masses that differ from their counterparts in Table HVl by more than 1.5cr. 



a (fm) 


P 




Mi 1 ^) 


Xco(l J Po) 


Xd(m) 




w 0.18 


6.503 


0.120 


4.213(1) 


★ 4.111(9) 


★ 4.210(9) 


4.272(15) 




6.485 


If 


4.227(7) 


4.105(8) 


★ 4.200(7) 


4.286(10) 




6.467 


If 


4.223(12) 


4.127(7) 


4.227(8) 


4.278(15) 




6.458 


If 


★ 4.253(9) 


4.128(9) 


★ 4.222(9) 


4.310(11) 


w 0.15 


6.600 


0.122 


4.492(7) 


4.344(6) 


4.458(7) 


4.537(11) 




6.586 




★ 4.493(7) 


4.349(7) 


4.462(7) 


4.536(13) 




6.572 


ff 


★ 4.516(9) 


4.375(9) 


4.488(9) 


4.574(10) 




6.566 


)) 


4.548(10) 


4.405(6) 


4.526(7) 


4.614(10) 


w 0.09 


7.11 


0.127 


5.199(11) 


5.030(8) 


5.170(12) 


5.257(12) 




7.09 


If 


★ 5.198(13) 


★ 5.034(11) 


★ 5.168(13) 


5.257(14) 




7.08 


ff 


5.178(15) 


★ 5.047(8) 


5.167(13) 


5.232(18) 



TABLE VI: Rest masses of the bottomonium states r/b, T, ht, Xbo> an d Xbi calculated with relativistic operators. All masses 
in units of n = 0.318 fm. The star denotes masses that differ from their counterparts in Table [VTT1 by more than 1.5cr. 



a (fm) 


P 


Kb 


r) b (l l S ) 




T(l a Si) 


T(2 a 5i) 


^(l'Pi) 


Xw(l"fl)) 


Xm(1 j A) 


» 0.15 


6.586 
6.572 


0.076 


7.3776(8) 
7.4061(9) 


8.202(5) 
8.241(63) 


7.4100(9) 
7.4386(9) 


8.209(5) 
8.248(7) 


8.269(120) 
8.321(13) 


★ 8.162(35) 
8.292(11) 


★ 8.147(40) 
8.318(11) 


» 0.12 


6.81 
6.79 
6.76, a 


0.086 


8.0690(10) 
8.0563(17) 
7.9815(9) 


8.933(12) 
8.910(14) 
8.870(11) 


8.1299(12) 
8.1167(19) 
8.0426(13) 


8.957(12) 
8.929(14) 
8.890(10) 


8.919(15) 
8.902(19) 
8.860(20) 


8.855(13) 
8.850(21) 
8.796(13) 


8.898(13) 
8.874(38) 
8.839(14) 


» 0.09 


7.11 
7.09 
7.08 


0.0923 

If 
ff 


10.2040(15) 
10.1861(11) 
10.1795(26) 


11.130(78) 
11.142(20) 
11.112(30) 


10.2627(19) 
10.2468(15) 
10.2397(33) 


11.160(32) 
11.161(18) 
11.137(56) 


11.050(14) 

★ 11.056(19) 

★ 11.066(19) 


★ 11.006(13) 
10.992(14) 

★ 11.017(12) 


★ 11.049(10) 
11.034(13) 

★ 11.048(14) 



TABLE VII: Rest masses of the charmonium states hb, Xbo, Xbi, an d Xbi calculated with nonrelativistic operators. All masses 
in units of ri = 0.318 fm. The star denotes masses that differ from their counterparts in Table ["VTl by more than 1.5a". 



a (fm) 


P 


K b 


hbil'Pi) 


Xbo(l'Po) 




Xf>2(l J P 2 ) 


» 0.15 


6.600 


0.076 


8.254(8) 


8.220(8) 


8.243(8) 


8.274(9) 




6.586 


If 


8.252(10) 


★ 8.216(10) 


★ 8.242(10) 


8.276(10) 




6.572 


If 


8.321(11) 


8.288(10) 


8.312(11) 


8.341(11) 




6.566 


If 


8.369(9) 


8.335(9) 


8.359(9) 


8.391(10) 


« 0.09 


7.11 


0.0923 


11.046(9) 


★ 10.981(10) 


★ 11.022(11) 


11.077(8) 




7.09 


ff 


★ 11.020(10) 


10.973(10) 


11.006(12) 


11.040(10) 




7.08 


ff 


★ 11.014(10) 


★ 10.964(11) 


★ 11.008(10) 


11.045(9) 
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A. Relativistic vs. nonrelativistic operators 

The statistical quality of our data can be judged from 
Fig. [1] which shows examples of typical two-point func- 
tions for the IS* pseudoscalar and vector states and their 
corresponding effective masses, calculated with relativis- 
tic operators. The data are from the coarse ensemble 
with ami/am s — 0.01/0.05. We have a clear signal and 
the effective masses have well-established plateaus. As 
already mentioned, for the IP states we used both rela- 
tivistic and nonrelativistic types of operators. Figure [2] 
compares the effective masses of the h c (lP) and hb(lP) 
states, calculated with both types of operators. The data 
show that the effective masses obtained with nonrelativis- 
tic operators plateau at an earlier i m i n . Despite the fact 
that the statistics in the nonrelativistic case are, in this 
example, three times lower than in the relativistic case, 
the errors on the fitted h c and hb masses are smaller than 
the ones calculated with relativistic operators. This find- 
ing holds for all of the IP states studied here. Our statis- 
tics with the nonrelativistic operators are 2-3 times lower 
than with the relativistic ones (except for the medium- 
coarse case where they are 6 times lower), yet the errors 
on the masses are up to 50% smaller, and in some cases 



smaller still — see Tables IIVHVIII for numerical compar- 
isons. The nonrelativistic operators couple much more 
weakly to the excited states and, thus, yield effective 
mass plateaus of better quality and fitted masses with 
smaller errors. All of our results for quarkonium masses 
are listed in Tables IIVHVIII with statistical errors calcu- 
lated with the bootstrap method and symmetrized. 

The central values of the IP states calculated with rel- 
ativistic and norelativistic operators occasionally differ 
by more than 1.5 uncorrelated a. This difference arises 
more often than expected, especially once correlations 
are considered. In the tables, these cases are labeled 
with a star. To check whether this difference is due to 
statistics, in some cases we carried out simultaneous fits 
to both the relativistic and nonrelativistic correlators. 
The masses extracted this way turned out to be indis- 
tinguishable from the masses from nonrelativistic data 
alone. This was not surprising, because the data from 
relativistic sources had larger fluctuations than that from 
nonrelativistic sources. Thus, in our further analysis of 
the chiral extrapolation and a dependence, we use the 
nonrelativistic results for the IP states wherever they 
are available. 
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B. k tuning in quarkonium and heavy-strange 
mesons 

In Sec. Ill CI we argued that the best way to tune the 
hopping parameter k is to use the spin-averaged kinetic 
mass of heavy-strange hadrons. If instead one would tune 
to the kinetic mass of the (spin-averaged) quarkonium 
ground state, the resulting tuned k could be different 
at nonzero lattice spacing. To study this discrepancy 
we have computed the quarkonium IS* kinetic mass for 
a wide range of k on the medium-coarse ensemble with 
ami/am s = 0.0290/0.0484. Figure El shows the results 
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FIG. 1: Propagators (a) and effective masses (b) for the r/ c , 
J/^}, r)b, and T states, with delta-function sources and sinks, 
from the coarse ensemble with ami/am s = 0.01/0.05. 



FIG. 2: Comparison between effective masses for h c and hb 
calculated with relativistic and nonrelativistic operators, on 
the fine ensemble with ami/am s — 0.0124/0.031. 
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and also shows the physical aM(lS) and aMy. 3 From 
a polynomial fit to the data, we get k c w 0.122 for the 
charmed quark, which is the same value as the one we ob- 
tain from matching to D s . However, because the relevant 
discretization effects arc larger in bottomonium than in 
B mesons, the tuned values of the hopping parameter 
differ substantially: Kb rs 0.094 from T vs. Kb « 0.076 
from B s . 

When we tune to the D s , some uncertainty in k arises. 
We take the tuning error in k c to be 0.0015 and in Kb 
to be 0.006. Reference [|| finds uncertainties (statistical 
and fitting) in this range on the medium-coarse, coarse, 
and fine ensembles, and here we assume the same for 
the extra-coarse ensembles. We discuss in the next sub- 
section how to propagate these errors to our computed 
splittings. 



TABLE VIII: l^i-l^o hyperfine splittings in ri units as a 
function of the valence k calculated for the medium-coarse 
ensemble with ami/am s = 0.0290/0.0484. 



K 


r 1 [M(l 3 Si)-M(l 1 S'o)] 


0.070 


0.0247(9) 


0.075 


0.0299(11) 


0.080 


0.0369(12) 


0.085 


0.0442(13) 


0.090 


0.0531(14) 


0.095 


0.0631(15) 


0.100 


0.0749(17) 


0.105 


0.0885(18) 


0.110 


0.1029(23) 


0.115 


0.1244(27) 


0.120 


0.1499(33) 


0.125 


0.1836(40) 


0.130 


0.2335(49) 



Above we mentioned a small difference in tuning the 
clover coupling for the coarse ensembles. The value of 
the tadpole coefficient uq used in that analysis was de- 
termined from mean Landau gauge link whereas the co- 
efficient used in the others was determined from the pla- 
quette. This difference means that our bare quark mass, 
i.e., k, has a slightly different definition on the coarse 
ensembles. Discrepancies in mass splittings caused by 
this choice should be eliminated via the nonperturbative 
tuning. 



jt-tuning uncertainties 
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FIG. 3: Spin-averaged kinetic mass aMi as a function of 
K, over a wide range, on the medium-coarse ensemble with 
ami/am s = 0.0290/0.0484. With a polynomial fit to the 
data, we can read off j(aM Vc + 3aMj/^) and aMr, finding 
k c w 0.122 and K b w 0.094. 



3 When this tuning was carried out, the rib had not yet been ob- 
served by experiment. 



Tables IIWvTlI and most of the plots in Sec. |IV] show 
statistical errors only, because the foremost aim of this 
paper is to understand the pattern of discretization er- 
rors. A systematic error also arises from inaccuracies in 
tuning k c and Kb, and to study the continuum limit it is 
necessary to propagate this error to the mass splittings. 
We discuss here how we treat these uncertainties. 

Several pieces of evidence show that the spin-averaged 
splittings depend very little on k. These splittings vary 
little from charmonium to bottomonium [24j, a feature 
understood to be a consequence of both systems ly- 
ing between the confining and Coulombic part of the 
potential [J This feature is, in fact, reproduced 
in our lattice-QCD data. Moreover, earlier work in 
the quenched approximation [Efij and with n/ = 2 | 
show negligible k dependence for spin-averaged split- 
tings. Thus, we shall assume that the K-tuning error 
for these splittings can be neglected. 

For spin-dependent splittings, we compute the IS* hy- 
perfine splitting as a function of k, on the medium-coarse 
ensemble with ami/am s — 0.0290/0.0484, the same en- 
semble as in Fig. [3] The data are summarized in Ta- 
ble EECO We fit the data to the form 



HFS 



I/k-I/kc, (3.1) 
bo/fi 2 + b,/^ 3 + 6 2 / M 4 + 6 3 //i 5 + b 4 /v 6 , (3.2) 



for K cr = 0.145, which enforces the requirement that, at 
large heavy quark mass tuq = fj,/2a, the splitting goes 
as 1/rriQ. The fit gives x 2 /dof = 0.6/8. From the fit 
result we estimate that an error of 0.0015 in the deter- 
mination of k c results in a 6% error in the charmonium 
hyperfine splitting, and an error of 0.006 in the determi- 
nation of Kb, a 22% error in the bottomonium hyperfine 
splitting. We expect that these errors are characteristic 
of all splittings driven by the spin-spin and tensor terms 
in the quarkonium effective potential, since in the non- 
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FIG. 4: The (a) \P-\S and (b) 1 P\-\S splittings in charmonium. The fine ensemble data are in blue fancy squares, the coarse 
in green circles, the medium-coarse in orange diamonds and the extra-coarse in red squares. The chirally extrapolated values 
are given in the legend and plotted with filled symbols. 



relativistic treatment, they all stem from the same term 
in the heavy-quark effective Lagrangian. 

The spin-orbit splitting remains to be considered. In 
our data and in experiment, it decreases from charmo- 
nium to bottomium similarly to the hyperfine and tensor 
splittings. Therefore, we shall assume the same relative 
error from the uncertainty in tuning k. 

Below we also present results for the splittings between 
twice the spin-averaged mass of D s and D*, and of B s 
and B*, and the corresponding IS* quarkonium mass. To 
estimate their K-tuning errors we have calculated these 
spin-averaged masses for several values of k near k c and 
Kb on the coarse ensemble with ami/am s = 0.01/0.05. 
These direct measurements allow us to propagate the k- 
tuning errors from the masses to the mass splittings. We 
obtain an error of 1.3% for charm and 13% for bottom. 
We assume the same error for these splittings at other 
lattice spacings. 



IV. SPECTRUM RESULTS 

We now present plots of quarkonium mass splittings as 
a function of the square of the sea-quark pion mass. The 
splittings and their errors are calculated using the boot- 
strap method. In most cases, we expect the dependence 
on the sea-quark mass to be mild, so we perform on our 
results a chiral extrapolation linear in M 2 down to the 
physical pion. The extrapolated values are denoted in 
each plot with filled symbols. The error bars come from 
symmetrizing the la (68%) interval of the bootstrap dis- 
tribution. 

Where possible, we compare our results to experi- 



mental measurements. As a rule we take the aver- 
age values from the compilation of the Particle Data 
Group The exception is the mass of the rjb(lS) 

meson, which has only recently been observed. We take 
M Vb — 9390.9 ± 2.8 MeV, based on our av erag e of two 
measurements by the BaBar Collaboration [48l [40(| and 
one by the CLEO Collaboration [50j]. 

In examining the results, we are interested in see- 
ing how well we can understand discretization errors via 
the nonrelativistic description of Eqs. (|2.13p - (|2.15p . We 
therefore carry out separate chiral extrapolations at each 
lattice spacing, and discuss whether the a dependence, 
and any deviations from experiment, make sense. 

From the effective Lagrangian discussion, we expect 
different discretization errors to affect spin-averaged and 
spin-dependent splittings. Errors in the spin-averaged 
splittings stem from the Darwin (D ■ E) term and the 
two p A terms. Errors in the spin-dependent splittings 
stem from the chromomagnetic (icr ■ B) and spin-orbit 
(icr ■ D x E) terms. Moreover, from the general structure 
of potentials arising from QCD [25|, [26j, we learn that 
icr B predominantly affects M (uSufs) an d M (nPt cnsor ) , 
while icr ■ D x E affects M(nP sp i n _ or bit)- 



A. Spin-averaged splittings 

Let us start with IP-IS and 1^-15 splittings, plot- 
ted in Figs. [4] and \5\ vs. (riM^) 2 . In the nonrelativistic 
picture, they arise predominantly at order v 2 via the ki- 
netic energy, which our tuning of k should normalize cor- 
rectly. The spin-dependent terms in C^q [cf. Eq. (|2. 15[) ] 
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open-charm threshold. The experimental point in (b) is not the spin-averaged splitting, but the T(25')-1S' mass difference, 
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do not contribute to spin averages (IS, IP) or to a spin 
singlet (1 Pi). Discretization errors remain, however, at 
order v 4 via the mismatches in Eqs. (|2.28|K(|2.30p . We 
assess these results using the error estimates in Ref . , 
which account for both the a dependence and the relative 
v 2 suppression. 

Our results for charmonium are shown in Fig. [4j Our 
results for both splittings approach the continuum phys- 
ical point as the lattice spacing decreases, and the size 
of the discretization effects is about what one expects: 



5-6% from mg ^ mi and 3-6% from 7714 7^ 7712 [5lj . 

The IP- 15 and \ X P\-\S splittings in bottomonium are 
given in Fig. [5] These splittings agree acceptably with ex- 
periment, given the estimated discretization errors, 2-3% 
from rriE 7^ f»2 and 2-5% from 777,4 ^ 7712 [51j . We can- 
not compare the /i(,(l 1 Pi) mass with experiment, because 
that state has not been observed [24[ , but our results for 
the l : Pi level agree very well with the l 3 Pj average. 

Next let us examine the 25-15 splitting. We fit a cor- 
relator matrix constructed from two interpolating opera- 
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FIG. 8: Splittings in bottomonium between the individual 25 states and the IS level. 



tors, local and smeared, to three or more states (i.e., two 
or more excited states). The error we assign to the mass 
determination estimates the uncertainties in our method. 
The results for charmonium as a function of (riM^) 2 are 
shown in Fig. [6^,. The lattice data appear to lie signifi- 
cantly above the experimental value at the smaller lattice 
spacings. The individual 25 levels show the same trends 
we observe in the spin-averaged level. In Fig. [7] we plot 
separately the ?7 C (25)-15' and V^S^-IS*. We see that 
both r] c (2S) and tp(2S) are responsible for the behavior 
seen in Fig. [6^, the latter especially so. The results for 
bottomonium (Fig. \Ejp) are more satisfactory. 



We suggest two possible reasons for the behavior of the 
charmonium 2S'-15 splitting results. First, the 25* are 
the only excited states in this study. Excited states are 
more difficult than ground states to determine accurately. 
With only two operators, our fits are less reliable, even 
though our fit model has at least three states. Second, 
the fit procedure does not take into account adequately 
the possible contribution of multiple open charm levels. 
For example, we have not used a two-body operator in the 
matrix correlator. With unphysically large quark masses, 
the open charm levels are unphysically high. As the sea 
quark mass is decreased, they come down. Moreover, the 
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box size of our lattices at the lightest sea quark mass 
is larger, which decreases the discrete level spacing of 
the would-be open-charm continuum. The dotted line in 
Fig. [5^ shows the location of the physical open charm 
threshold. It is dangerously close to the physical 25 lev- 
els, especially the ip(2S). Thus it is conceivable that 
nearby multiple open charm levels are being confused 
with the 2S and artificially raise its fitted mass. This 
explanation is consistent with the observed gradual rise 
of this level in the fine ensembles with decreasing light 
quark mass but not with the trends seen in the coarse 
and medium-coarse ensembles. 



For bottomonium in Fig. [SJa, the open bottom thresh- 
old is safely distant (off scale in this plot), so we do 
not expect a similar confusion in this channel. Figure [8j 
shows the individual 25 bottomonium levels separately. 
There is no comparison for the first excited pseudoscalar 
state i]b(2S)-lS, because the state has not yet been ob- 
served [24], although the extrapolated values appear to 
approach a consistent continuum limit. The first excited 
vector state splitting T(25)-15 is given in Fig. [5b. The 
chirally extrapolated values monotonically approach the 
experimental value and for the fine ensembles our split- 
ting agrees with the experiment. 
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FIG. 12: Spin-orbit splittings in IP levels, with M(lP sp i n _ or bit) defined in Eq. (|2.4p . for (a) charmonium and (b) bottomonium. 



B. Hyperfine splittings 

Now let us turn to the hyperfine structure. Our results 
for the hyperfine splitting in charmonium and bottomo- 
nium are presented in Figs. [9] and [TO] For the IS* levels 
and for 25* bottomonium, there is little dependence on 
the sea quark mass. To assess the approach to the con- 
tinuum limit one must bear in mind that the errors in 
Figs. l9l and fTOl are statistical only, and the systematic er- 
ror from K-tuning must also be taken into account. We 
thus take the values at the physical pion mass, apply the 



K-tuning error and plot these data vs. a 2 , as shown in 
Fig. [TTJ Both data sets are consistently linear in a 2 , so 
we carry out such an extrapolation. The extrapolated 
values in units of r\ are 0.187(12) for charmonium, with 
x 2 /dof = 1.9/2, and 0.087(20) for bottomonium, with 
X 2 /dof = 0.55/1. One can see, from comparing Fig. [TT1 
with Fig. [51 and 1 101 that the K-tuning uncertainties inher- 
ited from the heavy-strange kinetic mass are larger than 
the statistical uncertainties of the quarkonium rest-mass 
splittings. 

In physical units these extrapolated results are 
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Mj/^g) -M Vc(ls) = 116.0 ±7 A±™ MeV and M T(1S) - 
M^ (lg ) = 54.0 ± 12.4to o MeV, where the second er- 
ror comes from converting from r\ units to MeV. For 
charmonium the average of experimental measurements 
is 116.4 ± 1.2 MeV [241 ] . so our result is perfectly consis- 
tent. For bottomonium, the experimental measurements 
are 71.4±|;f ±2.7 MeV (H, 66.1^ ±2.0 MeV and 
68.5±6.6±2.0 MeV [13]; symmetrizing the error bars and 
taking a weighted average, we find M-wis) — M^ng) = 
69.4 ± 2.8 MeV. Our hyperfine splitting thus falls 1.2cr 
short. Note that with lattice NRQCD, the HPQCD Col- 
laboration finds M T{1S) - M rib{is) = 61 ± 4 ± 13 MeV 
(l9| , which agrees with the recent experimental measure- 



ments, yet also with our result. 

The errors on the final 15 hyperfine splittings quoted 
here encompass statistics (as amplified by extrapola- 
tions), k tuning, and r\. In addition, the coupling cb 
has been adjusted only at the tree-level, introducing an 
error of 0(a s a) that our continuum extrapolation would 
not eliminate. A preliminary result for the one-loop cor- 
rection to cb is available [56] , suggesting that a very small 
correction is needed beyond the tadpole improvement of 
Eq. (|2.26|) . when uo is set from the Landau link. 

The 25 hyperfine splittings for both charmonium and 
bottomonium are shown in Figs. [Hb and 110b . Unfortu- 
nately, these results are not very useful. Although the 



17 



2.0 



[Q 
S 

*7 1.8 



1.7 



-i 1 1 1 1 1 1 1 1 1 1 r 



s:i experiment 1 .749( 1 ) 
Kfine 1.770(12) 
o coarse 1.798(9) 
med. coarfle 1.876(17) 



-s #- 




_l I I I I I I I I I I l_ 



2.8 — 



S -2.6 



ICQ 

H. 2.4 
t, 



2.2 



1 1 1 1 1 1 1 1 1 1 
" S experiment 2.197(4) 


1 ' 


'_ Kfine 2.338(17) 




coarse 2.359(17) 




- <>med. coarse 2.559(37) 

T — — — 


- 


4— S Sr ~ - <l> — 


-# : 


1 i i i 1 1 i l 1 1 


1 . , : 



0.0 



0.5 

(r,Mj a 



1.0 



0.0 



0.5 

(r,Mj a 



1.0 



(a) (b) 

FIG. 15: Quarkonium-heavy-light splittings (a) 2M(Th) - M{TS) and (b) 2M(B~ S ) - M(1S). 




I g U I I I I I I I I I I I I I I I I I I I I l_ 

0.00 0.01 0.02 0.03 0.04 
a 2 (fm z ) 




J g Li i I i i i i I i i i i I i i i i I i i i i_ 

0.00 0.01 0.02 0.03 0.04 
a 2 (fm z ) 



(a) a linn (b) 

FIG. 16: Continuum extrapolations of (a) 2M(Th) - M(1S) and (b) 2M(B~ S ) - M(1S) 



charmonium splitting agrees, within large errors, with 
experiment, one should bear in mind the issue of thresh- 
old effects surrounding our determination of the ip(2S) 
mass, discussed above. The bottomonium splitting does 
not suffer from this problem, but the statistical and fit- 
ting errors are still too large to make a prediction of the 
as yet unobserved rjb(2S) mass. 



C. P-state splittings 

We now turn to splittings between the l 3 Pj levels, 
which stem from two contributions [25j]. As discussed 



above, one comes from exchanging a Coulomb gluon 
between a spin-orbit term, h^Hcr ■ (D x E)h^' in 
Eq. HUTS]), and the static potential, h^A 4 h^\ The 
other comes from exchanging a transverse gluon be- 
tween the chromomagnetic terms, hy r 'i(T ■ BhS + ^ and 
h(~Hcr-Bh(' . These two contributions can be separated 
by forming the combinations in Eqs. (|2.4p and (|2.5p [2(| . 

The spin-orbit splittings M(lP sp i n _ or t>it) are shown in 
Fig. [T2] for charmonium and bottomonium. They ex- 
hibit a small lattice-spacing dependence and agree well 
with experiment, indicating that the chromoelectric in- 
teractions and, hence, eg are adjusted accurately enough. 
The tensor splittings M(lP tensor ) are shown in Fig. [T5] 
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for charmonium and bottomonium. These chromomag- 
netic effects seem to approach the experimental value as 
a decreases. Since the tensor and the spin-spin poten- 
tial components both measure the effects of the chromo- 
magnetic interaction, we plot the ratio of the IP tensor 
splitting to the IS hyperfine splitting in Fig. 1141 The co- 
efficient mg 1 should drop out from the ratio and if there 
are no effects from higher-dimension operators, the ratio 
should be a constant which agrees with the continuum 
limit. If a higher-dimension operator has a significant 
contribution, then the ratio need not agree any better 
than the splittings themselves. The charmonium case, 
Fig. ITIk . seems to suggest that the higher-dimension op- 
erator matters, the bottomonium case, Fig. [14b. seems to 
suggest it does not. This outcome is plausible, because 
the v 2 suppression of the higher-dimension operator is 
10% in bottomonium, but only 30% in charmonium [cf. 
Eqs. (EHU) and 1(532)1 ] . 



D. Quarkonium vs. heavy-strange mesons 

Unlike other approaches to heavy quarks, lattice QCD 
is supposed to treat heavy-light mesons and quarkonium 
on the same footing. If we form the splitting 

2M(Dl) - M (IS) (4.1) 

the rest mass drops out, leaving a pure QCD quantity. 
Here M(D S ) denotes the spin average of D s and .D* 
masses. This mass difference is interesting from the point 
of view of the discretization effects, which should con- 
tribute less to the D s and B s than to the charmonium 
and bottomonium IS* states. We show this splitting (also 
for the bottom-quark sector) combining our quarkonium 
rest masses with the Fermilab-MILC heavy-strange rest 
masses @ in Fig. (TS] The correlation in the error is 
treated correctly with the bootstrap method, and, as 



elsewhere in this paper, the bootstrap errors are sym- 
metrized. Clearly, discretization effects are important at 
nonzero a. 

In Fig. 1161 we incorporate the K-tuning errors and 
show the a dependence of the above splittings. Carrying 
out an exptrapolation linear in a 2 , which is empirically 
suitable, wejind n[2M(p^) - M(TS)\ = 1.705 ± 0.021 
and n [2M (B s ) - M(1S)\_ = 2.19 ± 0.49; these corre- 
spond to 2M(D S ) - M (IS) = 1058 ± 13±g 4 MeV and 
2M(B~ S )-M(TS) = 1359 ±304^ MeV, with the uncer- 
tainty in 7*1 yielding the second error bar. The bottomo- 
nium extrapolation agrees with the experimental value, 
but the combined statistical and K-tuning errors are quite 
large. The charmonium extrapolation is la shy of the 
experimental value. Given the empirical nature of our 
continuum extrapolation, this is completely satisfactory. 



E. Summary of spectrum results 

To summarize our results, Fig. [17] shows the charmo- 
nium and the bottomonium spectra as splittings from the 
IS level and compares them to the experimental results. 
We have plotted the chirally extrapolated values at each 
lattice spacing and included statistical, K-tuning, and r\ 
uncertainties. Solid lines show the experimental values, 
where they are known, and dashed lines show estimates 
from potential models [54[ in other cases. 

For the splittings discussed above, Table HXl shows the 
continuum limit, taken via linear extrapolations in a 2 . 
One should bear in mind that the NRQCD-based the- 
ory of cutoff effects, explained in Sec. Ill C[ anticipates 
a less trivial lattice-spacing dependence. The linear-in- 
a 2 extrapolations are consistent with the data, which 
are not yet sufficient to resolve more complicated func- 
tional forms. In Table IIXI the second (asymmetric) er- 
ror bar comes from the conversion to MeV with r\ — 
0.318±g. O07 fm = 1.61118.035 GeV" 1 HHSI. 
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TABLE IX: Continuum extrapolations of splittings in charmonium and bottomonium in MeV. The first error comes from statis- 
tics and accumulated extrapolation systematics; the second comes from the uncertainty in scale setting with r% = 0.318lg qo? f m - 



Splitting 


Charmonium 


Bottomonium 




This work Experiment 


This work Experiment 



*Fi-lS 469±lliJ° 457.9 ±0.4 440 ± 17±o° 

25-TS 792 ±42±J 7 606 ± 1 599 ± 36±q 3 (580.3 ± 0.8) 

l'Si-l^o 116.0 ± 7.4±^ 6 116.4 ± 1.2 54.0 ± 12.4±^ 2 69.4 ± 2.8 

IP tensor 15.0 ± 2. 3l^ 3 16.25 ± 0.07 4.5 ± 2.2±° 1 5.25 ±0.13 

IP spin-orbit 43.3 ± 6.6±q 46.61 ± 0.09 16.9 ± 7.0±g 4 18.2 ± 0.2 

IS sQ-QQ 1058 ± 13lp 4 1084.8 ± 0.8 1359 ± 304lg 1 1363.3 ± 2.2 



a T(25)-lS instead of 25-15. 



The charmonium and bottomonium spectra by and 
large show good agreement with experiment. The char- 
monium hyperfine splitting agrees very well; the bot- 
tomonium splitting agrees at 1.2a. The tensor and spin- 
orbit splittings also agree well, for both systems. The 
■"Pi-lS and IP-IS spin-averaged splittings agree at 1.1- 
1.3cr for cc; the IP-IS at 0.6a for bb. As discussed above, 
the charmonium 2S states are too high, because our op- 
erator basis and statistics proved to be insufficient to 
disentangle the bound states from open-charm threshold 
effects. For bottomonium the 2S-1S splitting does not 
suffer from threshold effects and agrees well. When the 
r*i uncertainty is included, the splitting of quarkonium 
relative to the heavy-strange spectrum, 2M(D S ) — M(1S) 
and 2M(B S ) — M(1S), also agrees well with experiment. 

V. CONCLUSIONS 

Quarkonium properties offer an excellent test of lat- 
tice QCD, because they are relatively well-understood 
hadrons, via potential models and effective field theo- 
ries. This paper attempts a thorough study of the char- 
monium and bottomonium mass splittings, using lattice 
gauge fields with 2±1 flavors of sea quarks. By using the 
Fermilab method for heavy quarks, we are able to study 
both systems, as well as heavy-light hadrons, with the 
same basic theoretical tool. By using the MILC ensem- 
bles, we are able to study a wide range of lattice spacing, 
and a wide range of up and down sea-quark masses, down 
to 0.10m 5 . 

Our aim here has been to develop methods and to com- 
pare discretization effects against expectations that are 
gleaned from an effective theory analysis. An important 
technical finding for ground P states is that nonrelativis- 
tic operators are superior to relativistic operators in over- 
lap and, hence, statistics. 

Our calculations reproduce most features of the mass 
splittings, to the extent expected. This optmistic conclu- 
sion is marred somewhat, because we find that the errors 



from k tuning are significant for spin-dependent split- 
tings. Agreement with experiment is found only when 
these uncertainties, which stem from the heavy-strange 
kinetic mass, are taken in to account. In some other 
cases, such as leptonic decay constants for heavy-light 
mesons [111 ] , uncertainties in k also influence significantly 
the final error budget. 

In the continuation of this project, we hope to im- 
prove on the results presented here in several ways. 
First, the MILC ensembles now contain approximately 
four times as many configurations, and they extend to 
smaller lattice spacings, a w 0.06 fm and a « 0.045 fm. 
The finer lattice will bring charm into the region where 
Symanzik-motivated continuum extrapolations are justi- 
fied and should bring bottomonium discretization effects 
under 1%. To this end it may also prove worthwhile to 
incorporate the p 4 corrections of the improved Fermilab 
action 15111 . Higher statistics and twisted-boundary con- 
ditions [571 ] should improve the tuning of n and, thus, 
reduce errors from this source as well. 
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